1.1 Protein-protein interaction network description and loading

The first dataset consist in an interaction network table available on STRING in which each protein-protein interaction is annotated with different confidence scores. These scores should be interpreted as “indicators” reflecting the probability that the reported interaction is considered true given the existing evidences. So, for example, if the scores rank from 0 (lowest confidence) to 1 (highest confidence), a score of 0.5 would indicate that every second interaction might be a false positive.

#load E.Coli protein-protein interaction network
EcoPPI<-read.table("511145.protein.physical.links.full.v11.0.txt",header=TRUE)
head(EcoPPI)

##   protein1 protein2 homology experiments experiments_transferred database
## 1    b0002    b0003        0           0                     122        0
## 2    b0002    b0004        0           0                     122        0
## 3    b0002    b0008        0           0                       0        0
## 4    b0002    b0014        0           0                      92        0
## 5    b0002    b0026        0           0                       0        0
## 6    b0002    b0031        0           0                       0        0
##   database_transferred textmining textmining_transferred combined_score
## 1                    0          0                    477            521
## 2                    0          0                    539            577
## 3                    0          0                    128            128
## 4                    0          0                    152            197
## 5                    0          0                    209            209
## 6                    0          0                    540            540

cat(paste("The number of protein-protein interactions is", nrow(EcoPPI), 
      "and there are", ncol(EcoPPI)-2, "types of scores."))

## The number of protein-protein interactions is 368774 and there are 8 types of scores.

From the table above, is possible to notice that in most cases there are two types of scores: ‘normal’ and ‘transferred’. The difference between the two regards the fact that the latter is computed from data coming from other organism (and not the organism of interest) and then transferred via homology/orthology.

1.2 Gene expression compendium description and loading

The second dataset is the E.coli gene expression compendium, built using RNA-seq libraries from SRA in an automatic fashion.

#load E.coli gene expression compendium
load("Ecoli_compendium.RData")

It contains three different versions of gene expression data:

• ALL_FINAL_Nreads.

• ALL_FINAL_TPM.

• ZL2.

In all three cases, the rows of the matrices are genes (with the gene identifier) and the columns are the different conditions where gene expression is measured.

For the subsequent analysis, among the three, the ZL2 is chosen and merged with the information contained in the protein-protein interaction network.

Ecoli_comp <- as.data.frame(ZL2)
#merging PPI and ZL2 
Ecoli_comp$protein <- row.names(Ecoli_comp)
EcoPPI_GeneExpr <- merge(EcoPPI, Ecoli_comp, by.x="protein1", by.y="protein")

2.1 Confidence based filtering

As previously mentioned, in the PPI network, the presence of confidence scores makes possible to quantify the strength of the supporting evidence for each reported interaction.

In order to increase the strength of evidence and to avoid taking into consideration false-positive interactions, we filter the PPI network by imposing a score threshold both to the experiment score and the transferred experiment score.

Network selection

The experiment_transferred network is bigger and more connected (112328 interactions) with respect to the experiment one (4649 interactions).

Taking into account the computational cost, the confidence value distribution and the nature of the interactions, for the subsequent analysis we will focus on the network containing the experimentally verified interactions only.

2.3 Graph generation

After the filtering confidence procedure, it is possible to transform the network into a graph, which is an abstract object which describes a biological network with nodes and edges.

#graph creation
graph_eco_exp <-graph.data.frame(EcoPPI_EXP_filtered)

#extract connected components
connected_components <-components(graph_eco_exp)
cat(paste("The number of connected components is: ", length(unique(connected_components$membership))))

## The number of connected components is:  86

We can have a visual inspection of the graph by using:

• iGraph plot built-in function;

• Cytoscape

In both the two visualizations it is possible to notice the presence of a large number of interaction aggregated in the center of the graph.

2.4 Degree distribution

As mentioned above, a graph is made of two main elements, edges (links) and nodes, which can be used to analytically describe a network. In particular, we can define the degree of a node as the number of associated links to each single node and, starting from it we can calculate the degree distribution as the probability that a node has exactly k link. It is displayed as a density probability and the shape of the curves gives information about the TYPE of the network.

Based on the shape of the degree distribution we can define:

• RANDOM NETWORK in which each pair of nodes is connected with equal probability following a normal distribution (bell shaped).

• SCALE FREE NETWORK in which most of the nodes have few links (small probability) and only a small portion of nodes is connected to many other nodes (high probability). In that case we talk about power-law distribution.

Types of degree distribution

Most biological networks are scale-free, since they have only few important nodes connected to many other nodes (hubs) whereas the vast majority display few connections.

We can plot the degree distribution of our network and, since we are dealing with biological data, we expect that we will find a scale-free one.

degree_distr <- degree_distribution(graph_eco_exp, cumulative=T, mode="all")
plot(x=0:max(degree(graph_eco_exp)), y=degree_distr, pch=19, cex=1.2, col="#ADD8E6",
xlab="Degree (K)", ylab="Probability P(K)",
main = "Degree distribution")

Summarizing the findings we can say that the network’s degree distribution follows the power-law that characterizes scale-free networks and:

cat(paste("The network contains",vcount(graph_eco_exp),"proteins and", ecount(graph_eco_exp),"interactions."))

## The network contains 1188 proteins and 4649 interactions.

Hubs definition

As did by Haiyan et al. first of all we identify hubs as all the proteins that have the 20% highest number of neighbors.

#Identify hubs as proteins with 20% of the total connections
deg <-igraph::degree(graph_eco_exp)

deg_tau<-quantile(deg,0.9)

We can have a visual inspection of the number of proteins classified as hubs and non hubs

h <- V(graph_eco_exp)$name[deg >= deg_tau]
nh <- V(graph_eco_exp)$name[deg < deg_tau]

h_nh <- data.frame('classification' = c('hubs', 'non_hubs'),
                'number' = c(length(h), length(nh)))

p4 <- ggplot(h_nh, aes(classification, number, fill = classification)) +
  geom_bar(position = 'dodge', stat = 'identity') +
  geom_text(aes(x = classification, label = number),position = position_dodge(width = 0.9), vjust =-0.1 , size = 3, col = 'black') +
  scale_fill_manual(values = c("#DDA0DD","#FFF8DC"),
                    name = 'Classification',
                    labels = c('Hubs','Non hubs'))+
  labs(title="Hubs non-hubs classification",x="Type", y = "Frequency")+
  scale_x_discrete(label = c('Hubs', 'Non hubs'))+
  theme_light()


p4 + theme(plot.title = element_text(face='bold'))

As expected, the vast majority of the proteins are classified as non-hubs.

Betweenness definition

We can then define betweenness as the proteins that are in the top 20% in terms of betweenness centrality.

bet<-betweenness(graph_eco_exp,directed=FALSE)
#filter based on the quantile threshold
bet_tau<-quantile(bet,0.9)

Identification of the four categories of nodes in a network

Using the degree and the betweenness as thresholds, we can identify all the nodes categories we are interested in.

## Select the bottleneck-hubs (BH) 
BH <-V(graph_eco_exp)$name[bet>=bet_tau & deg>=deg_tau]

## Select the on-bottleneck hubs (NBH)
NBH <-V(graph_eco_exp)$name[bet<bet_tau & deg>=deg_tau]

## Select the non-hub (BNH)
BNH <-V(graph_eco_exp)$name[bet>=bet_tau & deg<deg_tau]

## Select the Non bottleneck non hubs (NBNH)
NBNH <-V(graph_eco_exp)$name[bet<bet_tau & deg<deg_tau]

cat("We obtained:\n",
    length(BH), "Bottleneck-hubs (BH)\n",
    length(NBH), "Non bottleneck-hubs (NBH)\n",
    length(BNH), "Bottleneck-Non hubs (BNH)\n",
    length(NBNH), "Non bottleneck-Non hubs (NBNH)")

## We obtained:
##  81 Bottleneck-hubs (BH)
##  48 Non bottleneck-hubs (NBH)
##  38 Bottleneck-Non hubs (BNH)
##  1021 Non bottleneck-Non hubs (NBNH)

Also in this case we can plot the result of the classification.

final_class <- data.frame('classification' = c('Bottleneck hubs', 'Non-bottleneck hubs', 'Bottleneck non-hubs', 'Non-bottleneck non-hubs'), 'number' = c(length(BH), length(NBH), length(BNH), length(NBNH)))

p5 <- ggplot(final_class, aes(classification, number, fill = classification)) +
  geom_bar(position = 'dodge', stat = 'identity') +
  geom_text(aes(x = classification, label = number),
            position = position_dodge(width = 0.9), 
            vjust = -0.1, size = 3, col = 'black') +
  scale_fill_manual(values = c("#D8BFD8", "#B0E0E6", "#008080", "#F5F5F5"),
                    name = 'Classification')+
  labs(title="Nodes classification",x="Type", y = "Frequency")+
  scale_x_discrete(label = c('BH','NBH','BNH','NBNH'))+
  theme_light()

p5 + theme(plot.title = element_text(face='bold'))

As confirmation we can observe that “non-bottleneck-non hubs” is the largest category.

4.1 Pearson Correlation Coefficient calculation

To calculate the correlation coefficient (PCC) distribution for all the interactors we iterate over a loop in order to obtain a correlation matrix of all the proteins interacting with BH[n]/NBH[n] vs themselves in all combinations. We store the distribution of Pearson for each member of each category and then we put all the pearson values together by plotting an histogram. In order to avoid to give more “weight” to hubs with a larger number of interactors we take care to specify the breaks equal to the fixed intervals set above. As a last step we store the counts and the middle point of each interval, obtaining a matrix populated by the values of all members of the group.

4.2 PCCs comparison

We are interested in the PCC distribution for both types of hubs because in their work, from which this analysis takes inspiration, Han et al. found out that the average PCCs of hubs follow clearly a bimodal distribution in the compendium of yeast.

When we talk about bimodal distribution we refer to a distribution with two peaks, or two local maximus, that are points in which data points stop increasing and start decreasing.

The presence of two peaks in data indicates the existence of two different groups, that, in the paper mentioned above, are identified as two distinct population of party and date hubs.

Between the two distribution above, the NBH seems to follow a bimodal distribution but alone the visual inspection is not sufficient to establish the bimodality and so we perform a statistical test on both the distribution in order to determine it.

In particular we use the Hartigans’ dip test available in the multimode Rpackage, which tests the number of modes in a given distribution. In this test for Xi ∼ F, i.i.d., the null hypothesis is that F is a unimodal distribution (i.e. mod0 = 1), and the test alternative is non-unimodal at least bimodal.

#dip test for BH
modetest(pearsonBH$Pearsons, mod0 = 1, method = "HH")

## 
##  Hartigan and Hartigan (1985) dip test
## 
## data:  pearsonBH$Pearsons
## Dip = 0.048208, p-value = 0.986
## alternative hypothesis: true number of modes is greater than 1

Since the test p-value is almost equal to ~1 we have no evidence to accept the alternative hypothesis of multimodality, so the BH distribution is unimodal.

We perform the test also for the NBH:

modetest(pearsonNBH$Pearsons, mod0 = 1, method = "HH")

## 
##  Hartigan and Hartigan (1985) dip test
## 
## data:  pearsonNBH$Pearsons
## Dip = 0.076482, p-value = 0.398
## alternative hypothesis: true number of modes is greater than 1

In that case the p-value is lower than before but since it is almost equal to ~0.4, also in that case we have no evidence to reject the null hypothesis and we can state that the distribution is unimodal.

4.3 Party and Date hubs classification

Even if our results are in contrast with literature findings, we can still assume the presence of party and date hubs in both bottleneck and non-bottleneck distributions.

In particular, following the same approach used by Han et al., we can classify as party hubs those with relatively high average PCC and date hubs with low PCC.

In our case, since there is no bimodal distribution suggesting us a natural boundary for partitioning date hubs from party hubs, by looking at the two distributions above, we can visually establish a boundary in correspondence of PCC = 0.2 to divide each of the two categories of hubs in:

• date hubs if they have PCC < 0.2

• party hubs if they have PCC >= 0.2

Using this threshold we expect that the majority of BH will be classified as date hubs, while in NBH we expect a more balanced situation due to the fact that the distribution is slightly different.

To proceed with the classification we can calculate the number of interactors classified as party or date hub in accordance with the chosen PCC threshold.

4.4 Hubs composition comparison between BH and NBH

From the classification we notice, as expected, that for BH the vast majority of the hubs are date hubs, whereas in the case of NBH we have a more balanced situation and an higher number of party hubs.

In interaction network, protein bottleneck connect different functional clusters and so it is likely that bottleneck hubs have an higher tendency to be date hubs rather than party hubs. (https://bmcbioinformatics.biomedcentral.com/articles/10.1186/1471-2105-6-39)

We can further explore the relationship between party and date hubs with bottleneck and non-bottleneck hubs, with the aim of discover more about the relationship between them.

#extract bottleneck hubs
BH_2_hub <- as.data.frame(table(BH_2$hub))
BH_2_hub$category <- 'Bottleneck'

#extract non-bottleneck hubs
NBH_2_hub <- as.data.frame(table(NBH_2$hub))
NBH_2_hub$category <- 'Non-bottleneck'

#create the merged df with party and date annotation
BH_NBH_DP <- rbind(BH_2_hub,NBH_2_hub)
colnames(BH_NBH_DP) <- c('Hubs','n','Category')

Plot the frequency to have a better visualization:

p8 <- ggplot(BH_NBH_DP, aes(Category, n, fill = Hubs)) +
  geom_bar(position="dodge", stat = 'identity') +
  geom_text(aes(x = Category, label = n),
            position = position_dodge(width = 0.9), 
            vjust = -0.2, size = 3, col = 'black') +
  scale_fill_manual(values = c("#F08080", "#66CDAA"), name = 'Type of hubs') +
  labs(title="Party and date hubs in BH and NBH",x="Category", y = "Frequency")+
  theme_light()

p8 + theme(plot.title = element_text(face='bold'))

From this plot seems that the distribution of the date and party hubs in the two category is slightly different, in particular for what regards date hubs in BH.

In order to asses if there is a statistically significant difference on what we are observing, we build a contingency table to be used in the calculation of a chi-squared test, which is commonly used to compare the distribution of a categorical variable in a sample or a group with the distribution in another one.

#construct the contingency table
contingency_table <- as.table(rbind(c(58,19), c(23, 29)))
dimnames(contingency_table) <- list(Hubs = c("Date", "Party"),
                    Categories = c("Bottleneck", "Non-bottleneck"))
contingency_table

##        Categories
## Hubs    Bottleneck Non-bottleneck
##   Date          58             19
##   Party         23             29

We can compute the chi-squared setting the null hypothesis as: no difference (relationship) exists between the categorical variables taken into consideration (in our case between the hub types and the BH/NBHs)

print(chi_test <- chisq.test(contingency_table))

## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  contingency_table
## X-squared = 11.548, df = 1, p-value = 0.0006783

The p.value of the test is significant (<0.05) so we can reject the null hypothesis and say that there is a statistically significant relationship between the party/date hubs and the corresponding class to which they are associated.

In particular, we can use the correlation plot of the residuals in order to better understand the associations:

#correlation plot
corrplot(chi_test$residuals, is.cor = FALSE, col = COL2("BrBG"))

In this plot positive residuals (in green) specify a positive associations between the row and the column variable, on the contrary in brown are shown negative residuals which are negative associations.

From the image above, it’s evident that there is a strong positive association between party hubs and the non-bottleneck category and, even if less marked, between date hubs and bottleneck hubs.

5.1 Heatmaps

pheatmap(allpearsonsBH,
         cluster_cols = FALSE, 
         color = colorRampPalette(c("#FFFFFF", "#DDA0DD", "#800080"))(length(allpearsonsBH))
)

pheatmap(allpearsonsNBH,
         cluster_cols = FALSE, 
         color = colorRampPalette(c("#F5F5F5", "#5F9EA0", "#4682B4"))(length(allpearsonsNBH)))

5.2 UMAP and 3D-plot

Additionally, we can further investigate if there is some inner structure by apply a dimensionality reduction technique (UMAP) and visualize the result with a 3D plot.

umap_BH <- umap(allpearsonsBH,n_components = 3)

plot3d_BH <- plot3d(x=umap_BH[,1],y=umap_BH[,2],z=umap_BH[,3],
                     type='s',radius=.1,
                     col= "#DDA0DD",
                     box = FALSE, 
                     main= "3D-Bottleneck hubs plot", cex.main = 2,
                     xlab= "1stComponent", ylab="2ndComponent", zlab= "3rdComponent")
plot3d_BH
rglwidget()

umap_NBH <- umap(allpearsonsNBH,n_components = 3)


plot3d_nbh <- plot3d(x=umap_NBH[,1],y=umap_NBH[,2],z=umap_NBH[,3],
                     type='s',radius=.1,
                     col= "#5F9EA0",
                     box = FALSE, 
                     main= "3D-Non-bottleneck hubs plot", cex.main = 2,
                     xlab= "1stComponent", ylab="2ndComponent", zlab= "3rdComponent")
plot3d_nbh
rglwidget()

#transform them into df and assign a color for the category
allpearsonsBHdf <- as.data.frame (allpearsonsBH)
allpearsonsBHdf$hub <- "#DDA0DD"
allpearsonsNBHdf <- as.data.frame(allpearsonsNBH)
allpearsonsNBHdf$hub <- "#5F9EA0"
allPearsons <- rbind(allpearsonsBHdf,allpearsonsNBHdf)

umap_ALL <- umap(allPearsons,n_components = 3)

plot3d_all <- plot3d(x=umap_ALL[,1],y=umap_ALL[,2],z=umap_ALL[,3],
                     type='s',radius=.1,
                     col= allPearsons$hub,
                     box = FALSE, 
                     main= "3D-Non-bottleneck hubs plot",
                     xlab= "1stComponent", ylab="2ndComponent", zlab= "3rdComponent")
legend3d("topright", legend = c("BH", "NBH"),
      col =  c("#DDA0DD", "#5F9EA0"), pch = 16, cex=1, inset=c(0.02))

plot3d_all
rglwidget()